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I. INTRODUCTION 

The state of a classical composite system is described in the phase space as a product 
of its individual constituents separate states. On the other hand, the state of a composite 
quantum system is expressed in the Hilbert space as a superposition of tensor products of its 
individual subsystems states. This means that the state of the quantum composite system is 
not necessarily expressible as a product of the individual quantum subsystems states. This 
oeculiar property of quantum systems is called Entanglement, which has no classical analog 
l|. The phenomenon of entanglement was first introduced by Schrodinger [2] who called 
it " Verschrankung" and stated: "For an entangled state the best possible knowledge of the 
whole does not include the best possible knowlege of its parts" . Quantum entanglement is 
a nonlocal correlation between two (or more) quantum systems such that the description 
of their states has to be done with reference to each other even if they are spatially well 
separated. In the early days of the quantum theory the notion of entanglement was first 
noted and introduced by Einstein et al. [3| as a paradox in the formalism of the quantum 
theory. Einestien, Podolsky and Rosen in their famous EPR paper proposed a thought 
experiment to demonstrate that the quantum theory is not a complete physical theory as it 
lacks the elements of reality needed for such a theory. It needed about three decades before 
performing an experiment that invalidated the EPR argument and guaranteed victory to 
the quantum theory. The experiment was based on a set of inequalities derived by John 
Bell [4], which relate correlated measurements of two physical quantities that should be 



obeyed by any local theory He demonstrated that the outcomes in the case of quantum 
entangled states violate the Bell inequality. This result emphasizes that entanglement is a 
quantum mechanical property that can not be simulated using a classical formalism. It took 
about three decades for the theoretical results confirming the Bell inequality violation to be 
obtained |5|. 

Recently the interest in studying quantum entanglement was sparked by the development 
in the field of quantum computing which was initiated in the eighties by the pioneering work 
of Benioff, Bennett, Deutsch, Feynman and Landauer |6l-ll2l|. This interest gained a huge 
boost in 1994 after the distinguishable work of Peter shore where he developed a quantum 
computer algorithm for efficiently prime factorizing the composite integers 



where entanglement plays a major role are quantum teleportatio n [14 



13j. Other fields 

n 

15|, dense coding 



IS Il7j], quantum communication [18] and quantum cryptography [191 ] . 



or 



;he underlying 



28]. The basic 



Different physical systems have been proposed as reliable candidates 
technology of quantum computing and quantum information processing [20 
idea in each one of these systems is to define certain quantum degree of freedom to serve as 
a qubit, such as the charge, orbital, or spin angular momentum. This is usually followed by 
finding a controllable mechanism to form an entanglement between a two-qubit system in 
such a way to produce a fundamental quantum computing gate such as an exclusive Boolean 
XOR. In addition, we have to be able to coherently manipulate such an entangled state to 
provide an efficient computational process. Such coherent manipulation of entangled states 
has been observed in different systems such as isolated trapped ions [29] and superconducting 
junctions [30[. The coherent control of a two-electron spin state in a coupled quantum dot 
was achieved experimentally, in which the coupling mechanism is the Heisenberg exchange 
interaction between the electron spins 3 lM33[] . 

Particularly, the solid state systems have been in the focus of interest as they facilitate 
the fabrication of large integrated networks that would be able to implement realistic quan- 
tum computing algorithms on a large scale. On the other hand, the strong coupling between 
a solid state system and its complex environment makes it a significantly challenging mis- 
sion to achieve the high coherence control required to manipulate the system. Decoherence 
is considered as one of the main obstacles toward realizing an effective quantum comput- 
ing system 34j-l371|. The main effect of decoherence is to randomize the relative phases of 
the possible states of the isolated system as a result of coupling to the environment. By 



randomizing the relative phases, the system loses all quantum interference effects and its 
entanglement character and may end up behaving classically. 

The interacting Heisenberg spin systems, in one, two and three-dimensions represent very 
reliable model for constructing quantum computing schemes in different solid-state systems 



and a very rich model for studying the novel physics of localized spin systems |38l44l|. These 
spin systems can be experimentally realized, for instance, as a one- dimensional chain and 
lattices of coupled nano quantum dots. 

Multi-particle systems are of central interest in the field of quantum information, as a 
quantum computer is considered as many body system by itself. Understanding, quantifying 
and exploring entanglement dynamics may provide an answer for many questions regarding 
the behaviour of complex quantum systems [42], particularly, quantum phase transitions 
and critical behaviour as entanglement is considered to be the physical property responsible 



for the long-range quantum correlations accompanying these phenomena 43H47]. 



A. Entanglement measures 

One of the central challenges in the theory of quantum computing and quantum infor- 
mation and their applications is the preparation of entangled states and quantifying them. 
There has been an enormous number of approaches to tackle this problem both experimen- 
tally and theoretically. The term entanglement measure is used to describe any function 
that can be used to quantify entanglement. Unfortunately we are still far from a complete 
theory which can quantify entanglement of a general multipartite system in pure or mixed 
state 48rol|]. There are very limited number of cases where we have successful entanglement 
measures. Two of these cases are (i) A bipartite system in a pure state and (ii) A bipartite 
system of two spin 1/2 in a mixed state. Since these two particular cases are of special 
interest for our studies in spin systems we discuss them in more details in the coming two 
subsections. 

To quantify entanglement, i.e. to find out how much entanglement is contained in a quan- 



tum state, Vedral et al. introduced the axiomatic approach to quantify entanglement 52 ] . 
They introduced the basic axioms that are necessary for an entanglement measure to sat- 
isfy. Before introducing these axioms one should first discuss the most common operations 



that can be performed on quantum systems and affect their entanglement [53] . (i) Local 



operation: is an operation that when applied to a quantum system consisting of two subsys- 
tems, each subsystem evolves independently. Therefore, possibly pre-existing correlations, 
whether classical or quantum, will not be affected. Hence, their entanglement also will not 
be affected, (ii) Global operation: is an operation that when applied to a quantum system 
consisting of two subsystems, the subsystems will evolve while interacting with each other. 
Therefore, correlations, both classical and quantum, may change under the effect of such op- 
erations. Hence, their entanglement will be affected by this operation, (iii) Local Operations 
with Classical Communications (LOCC): it is a special kind of global operations, at which 
the subsystems evolve independently with classical communications allowed between them. 
Information about local operations can be shared using the classical communications and 
then further local operations may be done according to the shared information. Therefore, 
classical, but not quantum, correlations may be changed by LOCC. It is reasonable then to 
require that an entanglement measure should not increase under LOCC. 

Now let us briefly introduce a list of the most commonly accepted axioms that a function 
E must obey to be considered an entanglement measure: (1) E is a mapping from the 
density matrix of a system to a positive real number p — > E(p) G R; (2) E does not increase 
under LOCC only [53]; (3) E is invariant under local unitary transformations; (4) For a 



pure state p = \ip) (?/>|, E reduces to the entropy of entanglement [54j (will be discussed in 



more details latter); (5) E = iff the states are separable 55j ; (6) E takes its maximum 



for maximally entangled states (Normalization); (7) E is continous 56J; (8) E should be a 



convex function, i.e. it cannot be obtained by mixing states [pi ^^,E(wiPi) < WiE(pi)]; (9) 
E is Additive, i.e. given two pairs of entangled particles in the total state a = o\ <g> o~2 then 
we have [57| E(a) = E(a{) + E(a 2 ). 

1. Pure bipartite state 

A very effective approach to handle pure bipartite state is the Schmidt decomposition. 
For a pure state \ip) of a composite system consists of two subsystems A and B, with two 
orthonormal basis {\4> a ,i)} and {\4>b,i)} respectively, the Schmidt decomposition of the state 
\ip) is defined by 

|V) =x>i^m>i<^>, (1) 



where Aj are positive coefficients satisfying J2 Xf = 1 and are called Schmidt coefficients. 
Evaluating the reduced density operators 

PA/B = tr B/A {\lp)) = ^2\ 2 i\<f)(a/b),i)(<P(a/b),i\ , (2) 

i 

shows that the two operators have the same spectrum Aj, which means that the two subsys- 
tems will have many properties in common. If the state \ip) is a cross product of pure states, 
say \(p a ,k) and \4>b,k), of the subsystems A and B respectively then \ip) is a disentangled state 
and all Schmidt coefficients vanish except A& = 1 in that case. In general the coefficients Aj 
can be used to quantify the entanglement in the composite system. 

Entropy plays a major role in the classical and quantum information theory. Entropy is 
a measure of our uncertainty (lack of information) of the state of the system. The Shannon 
entropy quantifies the uncertainty associated with a classical distribution {Pi} and is defined 
as H{x) = — Y^ x Px l°g Px- The quantum analog of the Shannon entropy is the Von Neumann 
entropy where the classical probability distribution is replaced with the density operators. 
Considering the density operator p representing the state of a quantum system, the Von 
Neumann entropy is defined as S(p) = — tr(plogp) = — ^c^logc^ where {a,} are the 
eigenvalues of the matrix p . For a bipartite system, the Von Neumann entropy of the reduced 
density matrix Pa/b namely S{pa/b) = ~J2i^i^°S^i * s a measure of entanglement which 
indeed satisfies all the above axioms. Although the entanglement content in a multipartite 
system is difficult to quantify, the bipartite entanglement of the different constituents of the 
system can provide a good insight about the entanglement of the whole system. 

2. Mixed bipartite state 

When the composite physical system is in a mixed state, which is more common, different 
entanglement measures are needed. In contrary to the pure state, which has only quantum 
correlations, the mixed state contains both classical and quantum correlations which the 
entanglement measure should discriminate between them. Developing and entanglement 
measure for mixed multipartite systems is a difficult task as it is very hard to discrimi- 



nate the quantum and classical correlations in that case 5J, [55|. Nevertheless for bipartite 
systems different entanglement measures were introduced which can overcome the math- 
ematical difficulty, particularly for subsystems with only two degrees of freedom. Among 



the most common entanglement measures in this case are the relative entropy Er 57 , 
Entanglement of distillation Ed 58]; Negativity and Logarithmic negativity 59]. 

One of the most widely used measures is entanglement of formation, which was the first 



measure to appear in 1996 by C. Bennett et al. [60]. For a mixed state the entanglement 
of formation is defined as the minimum amount of entanglement needed to create the state. 
Any mixed state \ip) can be decomposed into a mixture of pure states \ipi) with different 
probabilities Pi, which is called "ensemble" . The entanglement of formation of a mixed state 
can be obtained by summing the entanglements of each pure state after multiplying each 
one by its probability piEQtpi)). The entanglement of each pure state is expressed as the 
entropy of entanglement of that state. Hence, for an ensemble T of pure states {p i: \ipi)} we 
have [61] : 

E F (e) = J2PiEv N (m) (3) 

i 

Since a mixed state can be decomposed into many different ensembles of pure states with 
different entanglements, the entanglement of formation is evaluated using what is called the 
most "economical" ensemble 60|], i.e.: 

e f ( p ) = mJ2p* e M\^)), (4) 

i 

where the infimum is taken over all possible ensembles. Ep is called a convex roof and the 
decomposition leading to this convex roof is called the optimal decomposition 62]. 

The minimum is selected because if there is a decomposition where the average entan- 
glement is zero then this state can be created locally without the need for entangled pure 
states and therefore E F = 62]. 

Performing minimization over all decompositions is a very difficult task because of the 
large number of terms involved 5l| . Nevertheless, It was shown that only limited number 
of terms is sufficient to preform the minimization. However, finding explicit formula that 
does not need pref orming the minimization would simplify the evaluation of Ep significantly. 
Bennett et al. 60j evaluated Ep for a mixture of Bell's states, which are completely entangled 



qubits. Hill and Wootters [63j provided a closed form of Ep as a function of the density 
matrix for two-level bipartite systems having only two non-zero eigenvalues in terms of the 
concurrence^ which was extended latter to the case of all two-level bipartite systems, i.e. 



wo qubits 641 ] . The entanglement of formation satisfies the previously discussed axioms 



54, 



56J. 



B. Entanglement and quantum phase transitions 

Quantum phase transition (QPT) in many body systems, in contrary to classical phase 
transition, takes place at zero temperature. QPT are driven by quantum fluctuations as 



a consequence of Heisenberg uncertainty principle 65|, |66(. Examples of quantum phase 
transition are Quantum Hall transitions, magnetic transitions of cuprates, superconductor- 



insulator transitions in two dimension and metal-insulator transitions 



66 



671 ] . Quantum 



phase transition is characterized by a singularity in the ground state energy of the system as 



a function of an external parameter or coupling constant A 68]. In addition, QPT is char- 
acterized by a diverging correlation length £ in the vicinity of the quantum critical point 
defined by the parameter value A c . The correlation length diverges as £ _1 ~ J\X — X c \ v , 
where J is an inverse length scale and v is a critical exponent. Quantum phase transitions 
in many body systems are accompanied by a significant change in the quantum correlations 
within the system. This led to a great interest in investigating the behaviour of quantum 
entanglement close the critical points of transitions, which may shed some light on the dif- 
ferent properties of the ground state wave function as it goes through the transition critical 
point. On the other hand, these investigations may clarify the role entanglement plays in 
quantum phase transitions and how it is related to the different properties of that transi- 
tion such as its order and controlling parameters. The increasing interest in studying the 
different properties of entangled states of complex systems motivated a huge amount of re- 
search in that area. One of the main consequences of these research is the consideration of 
entanglement as a physical resource which can be utilized to execute specific physical tasks 



in many body systems 60|, |69| . Osborne and Nielsen have argued that the physical property 
responsible for long range quantum correlations accompanying quantum phase transitions 
in complex systems is entanglement, and that it becomes maximum at the critical point 



70|]. The renormalization group calculations demonstrated that quantum phase transitions 



have a universal character independent of the dynamical properties of the system and is 



only affected by specific global properties such as the symmetry of the system [TJj . In order 
to test whether the entanglement would show the same universal properties in the similar 
systems, the pairwise entanglement was studied in the XY spin model and its special case 



of the Ising model 44( where it was shown that the entanglement reaches a maximum value 



at the critical point of the phase transition in the Ising system. Also the entanglement was 



proved to obey a scaling behaviour in the vicinity of transition critical point for a class of 



one dimensional magnetic system, XY model in a transverse magnetic field 72]. Further- 
more, the scaling properties of entanglement in XXZ and XY spin-1/2 chains near and 
at the transition critical point were investigated and the resemblance between the critical 
entanglement in spin system and entropy in conformal field theories was emphasized [731 ] . 
Several works have discussed the relation between entanglement and correlation functions 
and as a consequence the notion of localizable entanglement was introduced, which enabled 
the definition of entanglement correlation length bounded from above by entanglement as- 
sistance and from below by classical functions and diverges at a quantum phase transition 
74H77J. Quantum discord which measures the total amount of correlations in a quantum 

state and discerns it from the classical ones, first introduced by Olliver and Zurek 781, 

I I 

was used to study quantum phase transition in XY and XXZ spin systems [79|. It was 

demonstrated that while the quantum correlations increases close to the critical points, the 
classical correlations decreases, in XXZ model, and is monotonous, in Ising model, in the 
vicinity of the critical points. 

C. Dynamics of entanglement 

In addition to the interest in the static behavior of entanglement in many body systems, 
its dynamical behaviour has attracted great attention as well, where different aspects of 
this dynamics have been investigated recently. One of the most important aspects is the 
propagation of entanglement through a many body system starting from a specific part 
within the system. The speed of propagation of entanglement through the system depends 
on different conditions and parameters such as the initial set up of the system, impurities 
within the system, the coupling strength among the system constituents, and the external 



magnetic field 80 



In most treatments, the system is prepared in an initial state described by an initial 
Hamiltonian Hi, then its time evolution is studied under the effect of different parameters, 
internal and external, which causes creation, decay, vanishing or just transfer of entanglement 
through the system. In many cases the system is abruptly changed from its initial state to 
another one causing sudden change in entanglement as well. 

The creation of entanglement between different parts of a many body system rather than 



10 



the transfer of entanglement through it was also investigated. The creation of entanglement 
between the end spins of a spin-1/2 XY chain was studied |83j. A global time-dependent 
interaction between the nearest neighbour spins on the chain was applied to an initial sepa- 
rable state and the creation of entanglement between the end spins on the chain was tested. 
It was demonstrated that the amount of entanglement created dynamically was significantly 
larger than that created statically. 

As heat can be extracted from a many body solid state system and be used to create heat, 
it was shown that entanglement can be extracted from a many body system by means of 
external probes and be used in quantum information processing [84]. The idea is to scatter 
a pair of independent non-interacting particles simultaneously by an entangled many body 
solid state system (for instance solid state spin chain) or optical lattice with cold atoms where 
each incident particle interacts with a different entangled system particle. It was demon- 
strated that the entanglement was extracted from the many body system and transferred 
to the incident probes and the amount of entanglement between the probes pair is propor- 
tional to the entanglement within the many body system and vanishes for a disentangled 
system. Recently the time evolution of entanglement between an incident mobile particle 
and a static particle was investigated [85]. It was shown that the entanglement increases 
monotonically during the transient but then saturates to a steady state value. The results 
were general for any model of two particles, where it was demonstrated that the transient 
time depends only on the group velocity and the wave packet width for the incident quais- 
monochromatic particle and independent of the type and strength of the interaction. On the 
other hand, entanglement information extraction from spin-boson environment using non- 



interacting multi-qubit systems as a probe was considered [86]. The environment consists of 
a small number of quantum-coherent two-level fluctuators (TLFs) with a damping caused 
by independent bosonic baths. A special attention was devoted to the quantum correlations 
(entanglement) that build up in the probe as a result of the TLF-mediated interaction. 

The macroscopic dynamical evolution of spin systems was demonstrated in what is known 
as Quantum domino dynamics. In this phenomenon, a one dimensional spin-1/2 system with 
nearest neighbour interaction in an external magnetic field is irradiated by a weak resonant 
transverse field 



87 



]. It was shown that a wave of spin flip can be created through the 
chain by an initial single spin flip. This can be utilized as a signal amplification of spin 
flipping magnetization. 

11 



II. DYNAMICS OF ENTANGLEMENT IN ONE-DIMENSIONAL SPIN SYS- 
TEMS 

A. Effect of a time-dependent magnetic field on entanglement 

We consider a set of N localized spin-| particles coupled through exchange interaction J 
and subject to an external magnetic field of strength h. We investigate the dynamics of en- 
tanglement in the systems in pre sence of a time- dependent magnetic fields. The Hamiltonian 



for such a system is given by 89| 



T N N N 

i=l i=l i=l 

where J is the coupling constant, h(t) is the time-dependent external magnetic field, a a are 
the Pauli matrices (a — x, y, z), 7 is the degree of anisotropy and N is the number of sites. 
We can set J=l for convenience and use periodic boundary conditions. Next we transform 
the spin operators into fermionic operators. So that, the Hamiltonian assumes the following 
form 

N/2 N/2 

H = } u a p (t) [CpC p + ct p c- p ] + i5 p [Cpct p + c v c„ v \ + 2h(t) = ^ H p . (6) 

p=i p=i 

where, a p (t) = —2cos<f) p — 2h(t), 5 P = 2 , ysin(j) p and p =27rp/N . It is easy to show [H p , H q ] = 

, which means the space of H decomposes into noninteracting subspace, each of four 

dimensions. No matter what h(t) is, there will be no transitions among those subspaces. 

Using the following basis for the pth subspace: (|0 >;Cpct p \0 >;Cp|0 >;c^ p |0 >), we can 

explicitly get 



H P (t) 



-»(*) 


— i8p 








Sp 


—4costf)p — 2h(t) 














—2cos4>p 






\ 



(7) 



V -2cos0p / 

We only consider the systems which at time t=0 are in the thermal equilibrium at tem- 
perature T. Let p p (t) be the density matrix of the pt\i subspace, we have p p (0) = e~ l3Hp ^ \ 
where (3 = 1/kT and k is the Boltzmann's constant. Therefore, using Eq. [3, we can 
have p p (0). Let U p (t) be the time-evolution matrix in the pth subspace, namely(/i = 1): 
i jjjQ = U p (t)H p (t) , with the boundary condition U p (0) = I . Now, the Liouville equation 
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C(i.k1) 



FIG. 1: Nearest-neighbor concurrence C at zero temperature as a function of the initial magnetic 
field a for the step function case with final field b. 



of this system is 



i d -f=m), *«)]. 



(8) 



which can be decomposed into uncorrelated subspaces and solved exactly. Thus, in the pth 
subspace, the solution of Liouville equation is p p (t) = U p (t)pp(0)Up(ty . 

As a first step to investigate the dynamics of the entanglement we can take the magnetic 
field to be a s tep function then generalize it to other relevant functional forms such as an 



oscillating one 



89[. Figure ([I]) shows the results for nearest-neighbor concurrence C(i,i + 1) 



at temperature T = and 7 = 1 as a function of the initial magnetic field a for the step 
function case with final field b. For a < 1 region, the concurrence increases very fast 
near b = 1 and reaches a limit C(i,i + 1) ~ 0.125 when b — > 00 . It is surprising that 
the concurrence will not disappear when b increases with a < 1. This indicates that the 
concurrence will not disappear as the final external magnetic field increase at infinite time. 
It shows that this model is not in agreement with the obvious physical intuition, since we 
expect that increasing the external magnetic field will destroy the spin-spin correlations 
functions and make the concurrence vanishes. The concurrence approaches maximum 



13 



C(i,i + 1) ~ 0.258 at (a = 1.37,6 = 1.37), and decreases rapidly as a ^ b. This indicates 
that the fluctuation of the external magnetic field near the equilibrium state will rapidly 
destroy the entanglement. However, in the region where a > 2.0, the concurrence is close 
to zero when b < 1.0 and maximum close to 1. Moreover, it disappear in the limit of b — > oo. 



Now, let us examine the system size effect on the entanglement with three different 
external magnetic fields changing with time t 90] : 



hi{t) 



h n {t) 



hin{t) 



b + (a — b)e 



-Kt 



a — asin(Kt) 







t< 
t> 

t <0 
t >0 



t < 
t > 



(9) 



(10) 



(11) 



a — acos(Kt) 

where a, b and K are varying parameters. 

We have found that the entanglement fluctuates shortly after a disturbance by an external 
magnetic field when the system size is small. For larger system size, the entanglement reaches 
a stable state for a long time before it fluctuates. However, this fluctuation of entanglement 
disappears when the system size goes to infinity. We also show that in a periodic external 
magnetic field, the nearest neighbor entanglement displays a periodic structure with a period 
related to that of the magnetic field. For the exponential external magnetic field, by varying 
the constant K we have found that as time evolves, C(i, i + 1) oscillates but it does not 
reach its equilibrium value at t — > oo. This confirms the fact that the nonergodic behavior 
of the concurrence is a general behavior for slowly changing magnetic field. For the periodic 
magnetic field hu = a(l — sin[—Kt\) the nearest neighbor concurrence is at maximum at 
t = for values of a close to one, since the system exhibit a quantum phase transition at 
A c = J/h = 1, where in our calculations we fixed J = 1. Moreover for the two periodic 
sin[—Kt] and cos[—Kt] fields the nearest neighbor concurrence displays a periodic structure 
according to the periods of their respective magnetic fields 90| . For the periodic external 
magnetic field hjjj(t), we show in Figure (|2J) that the nearest neighbor concurrence C(i,i + 1) 
is zero at t — since the external magnetic field hjjj(t = 0) = and the spins aligns along the 
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100 200 300 400 500 



FIG. 2: The nearest neighbor concurrence C(i,i + 1) (upper panel) and the periodic external mag- 
netic field hni{t) = a{\ — cos[Kt\), see Eq. (14) in the text (lower panel) for K = 0.05 with 
different values of a as a function of time t. 

x-direction: the total wave function is factorisable. By increasing the external magnetic field 
we see the appearance of nearest neighbor concurrence but very small. This indicates that 
the concurrence can not be produced without background external magnetic field in the Ising 
system. However, as time evolves one can see the periodic structure of the nearest nei ghb or 
concurrence according to the periodic structure of the external magnetic field hm(t) [90]. 



B. Decoherence in one dimesional spin system 

Recently, there has been a special interest in solid state systems as they facilitate the 
fabrication of large integrated networks that would be able to implement realistic quantum 
computing algorithms on a large scale. On the other hand, the strong coupling between a 



15 



solid state system and its complex environment makes it a significantly challenging mission 
to achieve the high coherence control required to manipulate the system. Decoherence is 
considered as one of the main obstacles toward realizing an effective quantum computing 



system 34H37J . The main effect of decoherence is to randomize the relative phases of the 
possible states of the isolated system as a result of coupling to the environment. By ran- 
domizing the relative phases, the system loses all quantum interference effects and may end 
up behaving classically. 

As a system of special interest, there has been great efforts to study the mechanism of 
electron phase decoherence and determine the time scale for such process (the decoherence 



ti me), in solid state quantum dots both theoretically [36ll9ll-l94J and experimentally [3ll.l33, 



95l-l97l|. The main source of electron spin decoherence in a quantum dot is the inhomogeneous 
hyperfine coupling between the electron spin and the nuclear spins. 

In order to study the decoherence of a two state quantum system as a result of coupling 
to a spin bath, we examined the time evolution of a single spin coupled by exchange inter- 
action to an environment of interacting spin bath modelled by the XY-Hamiltonian. The 



Hamiltonian for such system is given by 98] 



N 1 _ N N 

h = ~p E Ji,i + i*M + i - -f- E %* - E fc°? > ( 12 ) 

i=l i=l i=l 

where J^+i is the exchange interaction between sites i and % + 1, hi is the strength of the 
external magnetic field on site i, a a are the Pauli matrices (a — x,y,z), 7 is the degree of 
anisotropy and N is the number of sites. We consider the centered spin on the l t h site as the 
single spin quantum system and the rest of the chain as its environment, where in this case 
/ = (N + l)/2. The single spin directly interacts with its nearest neighbor spins through 
exchange interaction Jj_i^ = J^i+\ = J' ■ We assume exchange interactions between spins in 
the environment are uniform, and simply set it as J = 1. The centered spin is considered 
as inhomogeneously coupled to all the spins in the environment by being directly coupled 
to its nearest neighbors and indirectly to all other spins in the chain through its nearest 
neighbors. 



By evaluating the spin correlator C(t) of the single spin,the jth site 98] 

C,(t)=p|(t,/3)-pj(0,/3), (13) 

we observed that the decay rate of the spin oscillations strongly depends on the relative 
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FIG. 3: The spin correlation function C(t) of centered spin for N=501, h=0.5, and j = 1.0 versus 
time t for different values of the coupling (a) J' < J; (b) J' > J at zero temperature. The decay 
profile for each case is shown in the inner panel. 



magnitude of the exchange coupling between the single spin and its nearest neighbor J' and 
coupling among the spins in the environment J. The decoherence time varies significantly 
based on the relative couplings magnitudes of J and J'. The decay rate law has a Gaussian 
profile when the two exchange couplings are of the same order J' ~ J but converts to 
exponential and then a power law as we move to the regimes of J' > J and J' < J as shown 
in Fig. [3j We also showed that the spin oscillations propagate from the single spin to the 
environmental spins with a certain speed as depicted in Fig. HI Moreover, the amount of 
saturated decoherence induced into the spin state depends on this relative magnitude and 
approaches maximum value for a relative magnitude of unity. Our results suggests that 
setting the interaction within the environment in such a way that its magnitude is much 
higher or lower than the interaction with the single spin may reduce the decay rate of the 
spin state. The reason behind this phenomenon could be that the variation in the coupling 
strength along the chain at one point (where the single spin exits) blocks the propagation 
of decoherence along the chain by reducing the entanglement among the spins within the 

This 



environment which reduces its decoherence effect on the single spin in return 9! 

result might be applicable in general to similar cases of a centered quantum system coupled 

inhomogeneously to an interacting environment with large degrees of freedom. 
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FIG. 4: 77ie spin correlation function C(t) vs. time t with J' = 1.0 (black) and J 1 = 1.5 (red,) 
for the centered spin (upper left) and the spin at site 256 (upper right), 261 (lower left), and 266 
(lower right) in the environment. N=501, h=0.5, and 7 =1.0. 

C. An exact treatment of the system with step time-dependent coupling and 
magnetic field 

The obvious demand in quantum computation for a controllable mechanism to couple the 
qubits, led to one of the most interesting proposals in that regard which is to introduce a 
time-dependent exchange interaction between the two valence spins on a doubled quantum 
dot system as the coupling mechanism 
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39j . The coupling can be pulsed over definite 
intervals resulting a swap gate which can be achieved by raising and lowering the potential 
barrier between the two dots through controllable gate voltage. The ground state of the 
two-coupled electrons is a spin singlet, which is a highly entangled spin state. 

There has been many studies focusing on the entanglement at zero and finite temperature 
for isotropic and aniso trop ic Heisenberg spin chains in presence and absence of an external 
magnetic field 45|, |99h105| . Particularly, the dynamics of thermal entanglement has been 
studied in an XY spin chain considering a constant nearest neighbor exchange interaction, 
in presence of a time varying magnetic field represented by a step, exponential and sinusoidal 



89 



9Q| 



functions of time which we have discussed above 

Recently, The dynamics of entanglement in a one dimensional Ising spin chain at zer o 
temperature was investigated numerically where the number of spins was seven at most 106 ] . 
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The generation and transportation of the entanglement through the chain under the effect 
of an external magnetic field and irradiated by a weak resonant field were studied. It was 
shown that the remote entanglement between the spins is generated and transported though 
only nearest neighbor coupling was considered. Latter the anisotropic XY model for a small 
number of s pins, with a time dependent nearest neighbor coupling at zero temperature was 



studied too [1071 ] . The time- dependent spin-spin coupling was represented by a dc part and 
a sinusoidal ac part. It was found that there is an entanglement resonance through the chain 
whenever the ac coupling frequency is matching the Zeeman splitting. 

Here we investigate the time evolution of quantum entanglement in an infinite one di- 
mensional XY spin chain system coupled through nearest neighbor interaction under the 
effect of a time varying magnetic field h(t) at zero and finite temperature. We consider a 
time-dependent nearest neighbor Heisenberg coupling J(t) between the spins on the chain. 
We discuss a general solution for the problem for any time dependence form of the coupling 
and magnetic field and present an exact solution for a particular case of practical interest, 
namely a step function form for both the coupling and the magnetic field. We focused on 
the dynamics of entanglement between any two spins in the chain and its asymptotic be- 
havior under the interplay of the time-dependent coupling and magnetic field. Moreover, 
we investigated the persistence of quantum effects specially close to critical points of the 
system as it evolves in time and as its temperature increases. 

The Hamiltonian for the XY model of a one dimensional lattice with N sites in a time- 
dependent external magnetic field h{t) with a time-dependent coupling J{t) between the 
nearest neighbor spins on the chain is given by 

T(t) N T(t) N N 

8=1 4=1 j=l 

where oVs are the Pauli matrices and 7 is the anisotropy parameter. 

Following the stan dard procedure to treat the Hamiltonian (|14|) . we transform the Hamil- 



tonian into the form 



108 ] 



AT/2 

H = J2H P , (15) 

p =i 
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with H p given by 



->l n _J_ n* 



U 



# P = a p (t)[c p c p + c T _ p c_p] + iJ(t)5 p [cl&_ p + c p c_ p ] + 2/i(t) 



(16) 



where a p (t) = —2J(t) cos(f) p — 2h(t) and S p = 27sin0 p . 

Writing the matrix representation of H p in the basis {|0) , c p c_ p |0) , c P |0) , c_ p |0)} we 
obtain 

/ 2/i(t) -iJ{t)8 p \ 



-ffn 



(17) 



i.J{t)5 p -4J(t)cos(j) p -2h(t) 

-2J(t)cos0 p 

\ -2J(t)cos0 p / 

Initially the system is assumed to be in a thermal equilibrium state and therefore its 
initial density matrix is given by 



Pp(0) = e 



-/3H P (0) 



(18) 



where /3 = 1/kT, k is Boltzmann constant and T is the temperature. 

Since the Hamiltonian is decomposable we can find the density matrix at any time t, 
p p (t), for the pth subspace by solving Liouville equation given by 



ip p {t) = [H p {t),p p {t)} 



(19) 



which gives 



Pp (t) = U p (t) Pp (0)Ul(t) . (20) 

where U p (t) is time evolution matrix which can be obtained by solving the equation 



iU p {t) = U p {t)H p {t) . 
Since H p is block diagonal U p should take the form 



(21) 



U p (t) 



(22) 



fuf t Uf 2 \ 
UZx uh o o 

f/ 3 p 3 

\ o o o u p u j 

Fortunately, eq. (|2T|) may have an exact solution for a time-dependent step function form 
for both exchange coupling and the magnetic field which we adopt in this work. Other 
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time-dependent function forms will be considered in a future work where other techniques 
can be applied. The coupling and magnetic field are represented respectively by 



J(t) = J + (Ji - J )6(t) 



(23) 



h(t) = h + (hi - h )9(t) , (24) 

where 9(t) is the usual mathematical step function. With this set up, the matrix elements 
of Up can be evaluated. The reduced density matrix of any two spins is evaluated in terms 
of the magnetization defined by 

N „ 1/iV 



»4£w4^ 



(25) 



j=\ p=l 

and the spin-spin correlation functions defined by 



QX / QX QX \ QiJ / QlJ Qy \ QZ / QZ QZ \ 



(26) 



Using the obtained density matrix elemen ts on e can evaluate the entanglement between any 
pair of spins using the Wootters method [109 1. 



1. Transverse Ising Model 

The completely anisotropic XY model, Ising model, is obtained by setting 7 = 1 in the 
Hamiltonian (IT4"|) . Defining a dimensionless coupling parameter A = J/h, the ground state 
of the Ising model is characterized by a quantum phase transition that takes place at A close 
to the critical value A c = 1 [44]. The order parameter is the magnetization (a x ) which differs 
from zero for A > A c and zero otherwise. The ground state of the system is paramagnetic 
when A — > where the spins get aligned in the magnetic field direction, the z direction. For 
the other extreme case when A — > 00 the ground state is ferromagnetic and the spins are all 
aligned in the x direction. 

We explored the dynamics of the nearest-neighbor concurrence C(i, i + 1) at zero temper- 
ature while the coupling parameter (and or the magnetic field) is a step function in time. We 
found that the concurrence C(i, i + 1) shows a nonergodic behavior. This behavior follows 
from the nonergodic properties of the magnetization and the spin-spin correlation functions 
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FIG. 5: C(i, i + l) as a function of A for h = ho = hi and J = Jo = 3\ at (a) kT = with any combination 
of J and h; (b) kT = 1 with h = hi = 0.25, 1,4; (c) kT = 1 with J = Ji = 0.25, 1,4; (d) fcT = 3 with 
ho = hi = 0.25,1,4. 



as reported by previous studies 



89, 


110 



llll ] . At higher temperatures the nonergodic be- 



havior of the system sustains but with reduced magnitude of the asymptotic concurrence 
(as t — y oo). 

We studied the behavior of the nearest neighbor concurrence C(i,i + 1) as a function of 



A for different values of J and h at different temperatures. As can be seen in Fig. 5(a), the 



behavior of C(i,i + 1) at zero temperature depends only on the ratio J/h (i.e. A) rather 
than their individual values. Studying entanglement at non-zero temperatures shows that 
the maximum value of C(i, i + 1) decreases as the temperature increases. Furthermore, 



22 



C(i,i + 1) shows a dependence on the individual values of J and h, not only their ratio as 



illustrated in Fig. 5(b) , 5(c) and 5(d) We investigated the dependence of the asymptotic 





FIG. 6: The asymptotic behavior of C(i, i + 1) as a function of (a) Jo and J\ with ho = hi = 1; (b) hg and 
h\ with Jq — Ji — 2, where ho, hi and Jq are in units of J\, kT — and 7 = 1. 



behavior (as t — > oo) of the nearest neighbor concurrence on the magnetic field and coupling 



parameters ho, hi, J and J\ at zero temperature. In Fig. 6(a) we present a 3-dimensional 



plot for the concurrence versus Jq and J\ where we set the magnetic field at ho = h\ = 1, 



while Fig. 6(b) shows the asymptotic behavior of the nearest neighbor concurrence as a 
function of h and hi, while fixing the coupling parameter at J — J\ — 2. In both cases 
as can be noticed the entanglement reaches its maximum value close to the critical value 
A c = 1. Also as can concluded the behaviour of the asymptotic concurrence is much more 
sensitive to the change in the magnetic field parameters compared to the coupling ones. 

There has been great interest in investigating the effect of temperature on the quantum 
entangleme nt and th e critical behavior of many body systems and particularly spin systems 
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44j, 



68 



112 



1131 ] . Osborne and Nielsen have studied the persistence of quantum effects 



in the thermal state of the transverse Ising model as temperature increases 44j. Here we 
investigate the persistence of quantum effects under both temperature and time evolution 
of the system in presence of the time-dependent coupling and magnetic field. Interestingly, 
the time evolution of entanglement shows a very similar profile to that manifested in the 
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FIG. 7: The asymptotic behavior of C(i, i + 1) as a function of (a) A and kT, in units of Ji, with 7 = 1, 
ho = hi = 1 and Jq — J\; (b) Ai and kT with 7 = 1, /j,q = hi = 1 and Jo = 1. 

static case, i.e. the system evolves in time preserving its quantum character in the vicinity 
of the critical point and kT = under the time varying coupling. Studying this behavior 
at different values of Jq and /io shows that the threshold temperature, at which C(i,i + 1) 
vanishes, increases as A increases. 

2. Partially Anisotropic XY Model 



We now turn to the partially anisotropic system where 7 = 0.5. First, we studied the 
time evolution of nearest neighbor concurrence for this model and it showed a non-ergodic 
behavior, similar to the isotropic case, which also follows from the non-ergodic behavior of 
the spin correlation functions and magnetization of the system. Nevertheless, the equilibrium 
time in this case is much longer than the isotropic. We have investigated the behaviour of 
the nearest neighbour concurrence C(i, i + 1) as a function of A for different values of J and 
h and at different temperatures as shown in Fig. [HJ We first studied the zero temperature 
case at different constant values of J and h. For this particular case, C(i, % + 1) depends only 
on the ratio of J and h, similar to the isotropic case, rather than their individual values as 



depicted in Fig. 8(a) Interestingly, the concurrence shows a complicated critical behavior 
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FIG. 8: The asymptotic behavior of C(i, i + 1) with 7 = 0.5 as a function of A when ho — hi and Jo = J\ 
at (a) kT = with any combination of constant J and h; (b) kT = 1 with ho = hi = 0.25, 1,4; (c) kT = 1 
with J = Ji = 0.25, 1, 4; (d) kT = 3 with h = hi = 0.25, 1, 4. 

in the vicinity of A = 1, where it reaches a maximum value first and immediately drops to a 
minimum (very small) value before raising again to its equilibrium value. The raising of the 
concurrence from zero as J increases, for A < 1, is expected as in that case part of the spins 
which were originally aligned in the z-direction change directions into the x and y-directions. 
The sudden drop of the concurrence in the vicinity of A = 1, where A is slightly larger than 
1, suggests that significant fluctuations is taking place and the effect of J x is dominating 
over both J y and h which aligns most of the spins into the x-direction leading to a reduced 



25 



entanglement value. Studying the thermal concurrence in Fig. 8(b), 8(c) and 8(d) we note 
that the asymptotic value of C(i, i + 1) is not affected as the temperature increases. However 
the critical behavior of the entanglement in the vicinity of A = 1 changes considerably as the 
temperature is raised and the other parameters are varied. The effect of higher temperature 



is shown in Fig. 8(d) where the critical behavior of the entanglement disappears completely 



at all values of h and J which confirms that the thermal excitations destroy the critical 
behavior due to suppression of quantum effects. The persistence of quantum effects as 





kT 0.5 







FIG. 9: The asymptotic behavior of C(i, i + 1) as a function of (a) A and kT, in units of J\, with 7 = 0.5, 
ho = hi = 1 and Jo = J\\ (b) Ai and kT with 7 = 0.5, ho = h\ = 1 and Jo = 1. 

temperature increases and time elapses in the partially anisotropic case is examined and 
presented in Fig. [9j As demonstrated, the concurrence shows the expected behavior as a 
function of A and decays as the temperature increases. As shown, the threshold temperature 
where the concurrence vanishes is determined by the value of A, it increases as A increases. 
The asymptotic behavior of the concurrence as a function Ai and kT is illustrated. Clearly 
the non-zero concurrence shows up at small values of kT and Ai. The concurrence has two 
peaks versus Ai but as the temperature increases, the second peak disappears. 
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3. Isotropic XY Model 

Now we consider the isotropic system where 7 = (i.e. J x = J y ). We started with the 
dynamics of nearest neighbor concurrence, we found that C(i,i + 1) takes a constant value 
that does not depend on the final value of the coupling J\ and magnetic field h\. This 
follows from the dependence of the spin correlation functions and the magnetization on the 
initial state only as well. This is because the initial coupling parameters J x and J y , which 
are equal, force the spins to be equally aligned into the x- and y-directions, apart from those 
in the z-direction, causing a finite concurrence. Increasing the coupling parameters strength 
would not change that distribution or the associated concurrence at constant magnetic field. 
The time evolution of nearest neighbor concurrence as a function of the time-dependent 



(a) 









FIG. 10: (a) C(i, i + 1) as a function of Ai and t, in units of J 1 , at kT = with 7 = 0, ho = hi = 1 and 
Jo = 5; (b) C(i, i + 1) as a function of Ao and t at kT = with 7 = 0, ho = h\ = 1 and J\ = 5. 

coupling is explored in Fig. [TUJ Clearly, C(i, i + 1) is independent of Ai. Studying C(i, i + 1) 
as a function of Ao and t with ho = hi = 1 at kT = for various values of J\, we noticed 
that the results are independent of J\. Again, as can be noticed when J < ho the magnetic 
field dominates and C(i,i + 1) vanishes. While for J > h , C(i,i + 1) has a finite value, 
as discussed above. Finally we explored the asymptotic behavior of the nearest neighbor 
and next nearest neighbor concurrence in the A-7 phase space of the one dimensional XY 
spin system under the effect of a time-dependent coupling J(t) as shown in Fig. [TTJ As 
one can notice, the non- vanishing concurrences appear in the vicinity of A = 1 or lower and 
vanishes for higher values. One interesting feature is that the maximum achievable nearest 
neighbor concurrence takes place at 7 = 1, i.e. in a completely anisotropic system, while the 
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FIG. 11: (a) The asymptotic behavior of C(i, i + 1) as a function of Ai and 7 with ho = hi = 1 and Jo = 1 
at kT = 0; (b) The asymptotic behavior of C(i, i + 2) as a function of Ai and 7 with ho = hi = 1 and Jo = 1 
at fcT = 0. 

maximum next nearest neighbor concurrence is achievable in a partially anisotropic system, 
where 7 ~ 0.3. 

D. Time evolution of the driven spin system in an external time-dependent mag- 
netic field 

We investigate the time evolution of quantum entanglement in a one dimensional XY spin 
chain system coupled through nearest neighbor interaction under the effect of an external 
magnetic field at zero and finite temperature. We consider both time-dependent nearest 
neighbor Heisenberg coupling J(t) between the spins on the chain and magnetic field h(t), 
where the function forms are exponential, periodic and hyperbolic in time. Particularly we 
focus on the concurrence as a measure of entanglement between any two adjacent spins on 
the chain and its dynamical behavior under the effect of the time-dependent coupling and 
magnetic fields. We apply both analytical and numerical approaches to tackle the problem. 
The Hamiltonian of the system is given by 



8=1 4=1 j=l 

where <7;'s are the Pauli matrices and 7 is the anisotropy parameter. 



(27) 
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1. Numerical and exact solutions 



Following the standard procedure to treat the Hamiltonian (127J1 1081 ]. we again obtain 
the Hamiltonian in the form 

N/2 

H = J2H p , (28) 

P =i 

with H p given by 

H p = a p {t) [clc p + cl p c^ p ] + iJ{t)5 p [cld p + c p c. p ] + 2h{t) , (29) 

where a p (t) = —2J(t) cos0 p — 2h(t) and 5 P = 27sin0 p . Initially the system is assumed to 
be in a thermal equilibrium state and therefore its initial density matrix is given by 

P -PH P (0) 

P P (0) = ~^, Z = Tr(e-^), (30) 

where /3 = 1/kT, k is Boltzmann constant and T is the temperature. 

Since the Hamiltonian is decomposable we can find the density matrix at any time t, 
p p (t), for the pth subspace by solving Liouville equation, in the Heisenberg representation 
by following the same steps we applied in Eqs. fT9lT22l 

To study the effect of a time-varying coupling parameter J(t) we consider the following 
forms 

J exp (t) = Jx + (Jo-Ji)e- Kt , (31) 

Jcos(t) = Jo- Jo cos (Kt) , (32) 

■J S in(t) = Jo ~ Jo sin (Kt) , (33) 

tanh ( K(t - |) J +1 . (34) 



Jtanhif) — Jo + 



2 

Note that Eq. (|2T|) . in the current case, gives two systems of coupled differential equations 
with variable coefficients. Such systems can only be solved numerically which we carried out 
in this work. Nevertheless, an exact solution of the system can be obtained using a general 
time-dependent coupling J(t) and a magnetic field in the following form: 

J{t) = X h[t) (35) 

where A is a constant. Using (1211) and (I3"5"j) we obtain the non- vanishing matrix elements 

«ii« u W«u« n W5 -is P |J(()j (36) 

«21 ^22 / I u 2 \ U 2 2 J \ iS p -4 cos 4> p -\ 
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and 



iu 33 = -2 cos (f) p J(t) u 33 , M44 = m 33 . 
These system of coupled differential equations can be exactly solved to yield 

u u = cos 2 6e- iX ' £ J ^ dt ' + sin 2 Qe~^ & J ^ dt ' , 

-1 sin cos [e-^ /o '(*')*' - e -*> /o '<*'>*' } , 



W12 



M22 = sin 



W21 = -W12 , 

-tti J* J(f JcJf + cog 2 ^g-iAa /„' J(f )df 
3 2icos<£ p / * J(t')dt' 



■U33 = m 44 = e 



where 



sin 6* 



cos6> 



\ 



52 + (2 cos0 p + I ) - (2 cos0p + 



2,/(J2 + (2cos0p + 



\ 



yj% + (2 cos P + f ) + (2 cos P + f ) 
2 v /^ + (2cos0 p + f) 



Where the angles <fi and 6 1 were found to be 

= (72 + l)vr, tan 20 



2 cos P + f ' 



where n = 0, ±1, ±2, . . ., therefore 



sin 20 



cos 29 



2 cos d>v + 



P^ A 



(37) 

(38) 

(39) 
(40) 
(41) 
(42) 

(43) 

(44) 



(45) 



(46) 



181 + (2 cos P + f ) y^ + (2 cos P + § ) 

As usual, upon obtaining the matrix elements of the evolution operators, either numerically 
or analytically, one can evaluate the matrix elements of the two-spins density operator with 
the help of the magnetization and the spin-spin co rrela tion funtions of the system and finally 
produce the concurrence using Wootters method 1091 ] . 



2. Constant magnetic field and time-varying coupling 

First we studied the dynamics of the nearest neighbor concurrence C(i,i + 1) for the 
completely anisotropic system, 7 = 1, when the coupling parameter is J exp and the mag- 
netic field is a constant using the numerical solution. As can be noticed in Fig. [T2J the 
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FIG. 12: C(i,i + 1) as a function of t with J = 0.5, Ji = 2,h = 1,N = 1000 at kT = and (a) 
J = J exp , K = 0.1 (b) J - J exp , X - 10 ; (c) J = J tan/l , if = 0.1 ; (d) J = J tanh , K = 10. 



asymptotic value of the concurrence depends on K in addition to the coupling parameter 
and magnetic field. The larger the transition constant is, the lower is the asymptotic value 
of the entanglement and the more rapid decay is. This result demonstrates the non-ergodic 
behavior of the system, where the asymptotic value of the entanglement is different from 
the one obtained under constant coupling J x . 

We have examined the effect of the system size N on the dynamics of the concurrence, 
as depicted in Fig. [13j We note that for all values of iV the concurrence reaches an approx- 
imately constant value but then starts oscillating after some critical time t c , that increases 
as iV increases, which means that the oscillation will disappear as we approach an infinite 
one- dimensional system. Such oscillations are caused by the spin-wave packet propagation 
90(. We next study the dynamics of the nearest neighbor concurrence when the coupling 
parameter is J cos with different values of K, i.e. different frequencies, which is shown in 
Fig. [TH We first note that C(i,i + 1) shows a pe riodic behavior with the same period of 
J(t). We have shown in our previous work 114j that for the considered system at zero 
temperature the concurrence depends only on the ratio J/h. When J rj h the concurrence 
has a maximum value. While when J » h or J << h the concurrence vanishes. In Fig. [TH 
one can see that when J = J max , C(i,i + 1) decreases because large values of J destroy 
the entanglement, while C(i,i + 1) reaches a maximum value when J = J = 0.5. As J(t) 
vanishes, C(i,i + 1) decreases because of the magnetic field domination. 
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FIG. 13: C(i, i + 1) as a function of t (units of J" 1 ) with J = J exp , J = 0.5, J x = 2, /i = 1, K = 1000 at 
fcT = and AT varies from 100 to 300. 
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FIG. 14: Dynamics of nearest neighbor concurrence with 7 = 1 for J cos where Jo = 0.5, ft, = 1 at kT = 
and (a) K = 0.1 ; (b) X = 0.5. 

5. Time- dependent magnetic field and coupling 

Here we use the exact solution to study the concurrence for four different forms of coupling 
parameter J exp , Jtanh, Jcos and J s « n when J(t) = \h(t) where A is a constant. We have com- 
pared the exact solution results with the numerical ones and they have shown coincidence. 
Figure [151 show the dynamics of C(i,i + 1) when h(t) = J(t) = J cos and h(t) = J(t) = J sin 
respectively, where Jo = 0.5 and K — 1. Interestingly, the concurrence in this case does 
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not show a periodic behavior as it did when h(t) = 1, i.e. for a constant magnetic field, in 
Fig. [TU In Fig. ITBT a) we study the behavior of the asymptotic value of C(i, i + 1) as a func- 
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FIG. 15: Dynamics of nearest neighbor concurrence with 7 = 1 at kT = with J = h = 0.5, K = 1 for 
(a) Jcos and h cos ; (bj J S m and h S i n . 

tion of A at different values of the parameters Jo, J\ and K where J(t) = Xh(t). Obviously, 
the asymptotic value of C(i, i + 1) depends only on the initial conditions not on the form or 
behavior of J(t) at t > 0. This result demonstrates the sensitivity of the concurrence evo- 
lution to its initial value. Testing the concurrence at non-zero temperatures demonstrates 
that it maintains the same profile but with reduced value with increasing temperature as 
can be concluded from Fig. [TBIf b). Also the critical value of A at which the concurrence 
vanishes decreases with increasing temperature as can be observed, which is expected as 
thermal fluctuations destroy the entanglement. Finally, in Fig. [T7] we study the partially 
anisotropic system, 7 = 0.5, and the isotropic system 7 = with J(0) = 1. We note that 
the behavior of C{i ,i + 1) in this case is similar to the case of constant coupling parameter 
studied previously [114j . We also note that the behavior depends only on the initial cou- 
pling J(0) and not on the form of J(t) where different forms of J(t)have been tested. It is 



interesting to notice that t 
of the previous works [89, 



he resu 



90. 



ts of Figs. [151 [16] and H7] confirm one of the main results 
114| namely that the dynamic behavior of the spin system, 
including entanglement, depends only on the parameter A = J/h not the individual values 
of h and J for any degree of anisotropy of the system. In these previous works both the 
coupling and magnetic field were considered time-independent, while in this work we have 
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FIG. 16: The behavior asymptotic value of C(i,i + 1) as a function of A with j = 1 at (a) kT — ; (b) 
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FIG. 17: The behavior asymptotic value of C(i,i + 1) as a function of A at kT = with (a) 7 = 0.5 ; (b) 

7 = 0. 

assumed J(t) = Xh(t) where h(t) can take any time-dependent form. This explains why 
the asymptotic value of the concurrence depends only on the initial value of the parameters 
regardless of their function form. Furthermore in the previous works, it was demonstrated 
that for finite temperatures the concurrence turns to be dependent not only on A but on the 
individual values of h and J, while according to Fig. HW b) even at finite temperatures the 
concurrence still depends only on A where J(t) = XJ(t) for any form of h(t). 



34 



III. DYNAMICS OF ENTANGLEMENT IN TWO-DIMENSIONAL SPIN SYS- 
TEMS 

A. An exact treatment of two-dimensional transverse Ising model in a triangular 
lattice 



Entanglement close to quantum phase transitions was originally analyzed by Osborne and 
Nielsen 44 1. and Osterloh et al. [72] for the Ising model in one dimension. We studied before 



a set of localized spins coupled through exchange interaction and subject to an external 



magnetic filed 46|, 



m 



89 



90] . We demonstrated for such a class of one-dimensional magnetic 



systems, that entanglement can be controlled and tuned by varying the anisotropy parameter 
in the Hamiltonian and by introducing impurities into the systems. In particular, under 
certain conditions, the entanglement is zero up to a critical point A c , where a quantum 
phase transition occurs, and is different from zero above A c [421 ] . 

In two and higher dimensions n early all calculations for spin systems were obtained 
by means of numerical simulations 1151 ] . The concurrence and localizable entanglement 
in two -dim e nsio nal quantum XY and XXZ models were considered using quantum Monte 
Carlo llq Ill7] . The results of these calculations were qualitatively similar to the one- 
dimensional case, but entanglement is much smaller in magnitude. Moreover, the maximum 
in the concurrence occurs at a position closer to the critical point than in the one- dimensional 
case [62]. 



118. 


119 



In this section, we introduce how to use the Trace Minimization Algorithm 
to carry out an exact calculation of entanglement in a 19-site two dimensional transverse 
Ising model. We classify the ground state properties according to its entanglement for 
certain class on two-dimensional magnetic systems and demonstrate that entanglement can 
be controlled and tuned by varying the parameter A = h/J in the Hamiltonian and by 
introducing impurities into the systems. We discuss the relationship of entanglement and 
quantum phase transition, and the effects of impurities on the entanglement. 



1. Trace minimization algorithm 

Diagonalizing a 2 19 by 2 19 Hamiltonian matrix and partially tracing its density matrix is 
a numerically difficult task. We propose to compute the entanglement of formation, first by 
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applying the trace minimization algorithm (Tracemin) |118l . Ill9| to obtain the eigenvalues 
and eigenvectors of the constructed Hamiltonian. Then, we use these eigenpairs and new 
techniques detailed in the appendix to build partially traced density matrix. 

The trace minimization algorithm was developed for computing a few of the smallest 
eigenvalues and the corresponding eigenvectors of the large sparse generalized eigenvalue 
problem 

AU = BUZ, (47) 

where matrices A,B G C nxn are Hermitian positive definite, U = [u\, ..., u p ] G <C nxp and X G 
MP xp is a diagonal matrix. The main idea of Tracemin is that minimizing Tr(X*AX), subject 
to the constraints X*BX = I, is equivalent to finding the eigenvectors U corresponding to 
the p smallest eigenvalues. This consequence of Courant-Fischer Theorem can be restated 



as 



iin_ Tr(X*AX) = Tr(U*AU) = ^ A;, (48) 



mm 

X*BX= 

i=l 



where I is the identity matrix. The following steps constitute a single iteration of the 
Tracemin algorithm: 

G = X* k BX k (compute G) 

G = VQV* (compute the spectral decomposition of G) 

H = Q*AQ (compute H, where Q = X fc V^- 1/2 ) 

H = WQW* (compute the spectral decomposition of H) 

X k = QW (now X* k AX k = A and X*BX k = I) 

X k+1 = X k - D ( D is determined s.t. Tr(X* k+1 AX k+1 ) < Tr(X* k AX k ) ). 
In order to find the optimal update D in the last step, we enforce the natural constraint 

X k BD = 0, and obtain 

' A BX ] ( D ) = ( M *\ (49) 

x;b o j \lJ \ o ) 

Considering the orthogonal projector P = BX k (XlB 2 X k )~ l XlB and letting D = (I — 
P)D, the linear system ( 1491) can be rewritten in the following form 

(/ _ p) A (I - p)D = (I- P)AX k . (50) 



Notice that the Conjugate Gradient method can be used to solve ( 1501) . since it can be shown 
that the residual and search directions r, p G Range(P)" 1 . Also, notice that the linear system 
(150]) need to be solved only to a fixed relative precision at every iteration of Tracemin. 
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FIG. 18: Diagonal elements of ^2 i a\ for N spins. 

A reduced density matrix, built from the ground state which is obtained by Tracemin, 
is usually constructed as follows: diagonalize the system Hamiltonian H(\), retrieve the 
ground state |\l/ > as a function of A = h/J, build the density matrix p = |\P >< \P|, and 
trace out contributions of all the other spins in density matrix to get reduced density matrix 
by p(i,j) = J2 P < u i(A)\ < v p (B)\p\uj(A) > \v p (B) >, where Ui(A) and v p (B) are bases 
of subspaces e(A) and e(B). That includes creating a 2 19 x 2 19 density matrix p followed 
by permutations of rows, columns and some basic arithmetic operations on the elements 
of p. Instead of operating on a huge matrix, we pick up only certain elements from \^f >, 
performing basic algebra to build a reduced density matrix directly. 

2. General forms of matrix representation of the Hamiltonian 



a* ® 



• / and J2i I ® • • • a l 



■L 



By studying the patterns of J^ <i ■> / (8) • • • of <g) 
one founds the following rules. 

a. Yli a i f or N spins The matrix is 2 N by 2^; it has only 2 N diagonal elements. 
Elements follow the rules shown in Fig. [18j 

If one stores these numbers in a vector, and initializes v = (N), then the new v is the 
concatenation of the original v and the original v with 2 subtracted from each of its elements. 
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We repeat this N times, i.e., 



v 
v-2 



(51) 




(52) 
(53) 



(54) 



*■ E<„>/ 



•07 



■afi 



, 10 / 1 , 

/ for N spins Since | ^ I exclude each 

1/ 10 



other, for matrix / (g) • • • of <g) • • • erf <g) • • • /, every row/column contains only one "1" , then 
the matrix owns 2 N "l"s and only "1" in it. If we know the position of "l"s, it turns out 
that we can set a 2 N by 1 array "col" to store the column position of "l"s corresponding to 
the 1st —¥ 2 Ar th rows. In fact, the non-zero elements can be located by the properties stated 
below. For clarity, let us number N spins in the reverse order as: N-l, N-2, . . . , 0, instead 
of 1, 2, . . . , N. The string of non-zero elements starts from the first row at: 1 + 2 % + 2 J ; with 
string length 2 J ; and number of such strings 2 J_1 . For example, Fig. [j~9l shows these rules 
for a scheme of / <E> <rf <g> o^ <8> I <8> I- 

Again, because of the exclusion, the positions of non-zero element "1" of I ® ■ ■ ■ af <g> 



• • • (Tj <8) • • • / are different from those of / <8> • • • erf ® • • • a^ <8> • • • /. So J^^„- ..^ / ® • • • of 



'<*j> 



o-" 



/ is a 2 by 2 matrix with only 1 and 0. 



After storing array "col", we repeat the algorithm for all the nearest pairs < i,j >, and 



concatenate "col"s to position matrix "c" of J^ 
section we apply these rules to our problem. 



<i,j> 



01 



■ ■ ■ o x a <8> • • • I- In the next 



3. Specialized matrix multiplication 



Using the diagonal elements array "v" of Yli °~i an d position matrix of non-zero elements 



" c " for E<i,i> /( 8)-- -of ®-- -o] 



/, we can generate matrix H, representing the Hamil- 
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FIG. 19: Scheme of matrix J <g> of <g) of <g) I <g) J 

tonian. However, we only need to compute the result of the matrix-vector multiplication 
H*Y in order to run Tracemin, which is the advantage of Tracemin, and consequently do 
not need to explicitly obtain H. Since matrix-vector multiplication is repeated many times 
throughout iterations, we propose an efficient implementation to speedup its computation 
specifically for Hamiltonian of Ising model (and XY by adding one term). 

For simplicity, first let Y in H*Y be a vector and J = h — 1 (in general Y is a tall matrix 
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FIG. 20: The illustration of H*Y. 



and J 7^ h ^ 1). Then 
H*Y 

(l 1 1 

1 1 1 



V 1 





( Y(l) \ 


/ 




Y(2) 




* 




+ 




\ Y(2 N ) j 


\ 



v(l) 



v(2) 



Y(2) 



1 V 

Y(c(l, 1)) + Y{c{l, 2)) + . . . + y(c(l, ^of pairs) 
Y(c(k, 1)) + F(c(A;, 2)) + . . . + Y(c(&, j^of pairs) 



v{2 N )J \Y(2 N )J 



\ 



( 



v(l)*Y(l) 
u(2)*Y(2) 



\ 



(,55) 



y Y (c(2 w , 1)) + y (c(2 N , 2)) + ... + Y(c(2 N , #of pairs) ) \ v(2 N ) * Y(2 N ) ) 

where p# stands for the number of pairs. 

When Y is a matrix, we can treat Y (2^ by p) column by column for ^<« 7 > I ® ■ • " a i ® 
■ ■ ■ crj (g) ■ • ■ /. Also, we can accelerate the computation by treating every row of Y as a vector 
and adding these vectors at once. Fig. [20] visualized the process. 

Notice that the result of the multiplication of the xth row of J2 <{ ,-> of crj (delineated by 
the red box above) and Y, is equivalent to the sum of rows of Y, whose row numbers are 
the column indecis of non-zero elements' of the xth row. Such that we transform a matrix 
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FIG. 21: Concurrence of center spin and its nearest neighbor as a function of A for both 7-site and 
19-site system. In the 7-site system, concurrence reaches maximum 0.15275 when A = 2.61. In the 
19-site system, concurrence reaches the maximum 0.0960 when A = 3.95. 

operation to straight forward summation & multiplication of numbers. 



4- Exact entanglement 

We examine the change of concurrence between the center spin and its nearest neighbor 
as a function of A = h/ J for both the 7-site and 19-site systems. In FigJJTJ the concurrence 
of the 7-site system reaches its maximum 0.152 75 w hen A = 2.61. In the 19-site system, 
the concurrence reaches 0.0960 when A = 3.95 120| . The maximum value of concurrence 
in the 19-site model, where each site interacts with six neighbors, is roughly 1/3 of the 
maximum concurrence in the one- dimensional transverse Ising mo del with s ize N=201 46] , 
where it has only two neighbors for each site. It is the monogamy 12ll . Il22j that limits the 
entanglement shared among the number of neighboring sites. This property also shown in 
the fact that fewer the number of neighbors of a pair the larger the entanglement among 
other nearest neighbors. Our numerical calculation shows that the maximum concurrence 
of next-nearest neighbor is less than 10~ 8 . It shows that the entanglement is short ranged, 
though global. 
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B. Time evolution of the spin system 



Decoherence is considered as one of the main obstacles toward realizing an effective quan- 



tum computing system 34] . The main effect of decoherence is to randomize th e re 



phases of the possible states of the considered system. Quantum error correction [123j and 



ative 



decoherence free subspace [35j, Il24| have been proposed to protect the quantum property 
during the computation process. Still offering a potentially ideal protection against environ- 
mentally induced decoherence is difficult. In NMR quantum computers, a series of ma gneti c 
pulses were applied to a selected nucleus of a molecule to implement quantum gates [125 ]. 
Moreover, a spin-pair entanglement is a reasonable measure for decoherence between the 
considered two-spin system and the environmental spins. The coupling between the sys- 
tem and its environment leads to decoherence in the system and vanishing of entanglement 
between the two spins. Evaluating the entanglement remaining in the considered system 
helps us to understand the behavior of the decoherence between the considered two spins 
and their environment 126 ] . 



The study of quantum entanglement in two-dimensional systems possesses a number 
of extra problems compared with systems of one dimension. The particular one is the 
lack of exact solutions. The existence of exact solutions has contribute d enormously to 



89 



108, 



1141 ]. Studies can be 



the understanding of the entanglement for ID systems [68 . 
carried out on interesting but complicated properties, applied to infinitely large system, and 
so forth use finite scaling method to eliminate the size effects, etc. Some approximation 
methods, li ke Density m atrix renormalization group (DMRG), are also only workable in one 



dimension 



125 



127H129J . So when we carry out this two-dimensional study, no methods can 



be inherited from previous researches. They heavily rely on numerical calculations, resulting 
in severe limitations on the system size and properties. For example, dynamics of the system 
is a computational-costing property. We have to think of a way to improve the effectiveness 
of computation in order to increase the size of research objects; mean while dig the physics 
in the observable systems. It may show the general physics, or tell us the direction of less 
resource-costing large scale calculations. 

To tackle down the problem, we introduce two calculation methods: step by step time- 
evolution matrix transformation and step by step projection. We compare them side by 
side, in short, besides the exactly same results, step by step projection method turned out 
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to be twenty times faster than the matrix transformation. 



1. The evolution operator 



According to quantum mechanics the transformation of \ipi(t )), the state vecto r at the 



initial instant to, into \i[)i(t)), the state vector at an arbitrary instant, is linear 1301 ] . There- 
fore there exists a linear operator U(t, t ) such that: 

\A(t))=U(t,t )\Mh))- (56) 

This is, by definition, the evolution operator of the system. Substituting Eq. (1561) into the 
Schorodinger equation, we obtain: 

ih^U(t,t )m )) = H(t)U(t,t )\i>(t )), (57) 

which means 

m^U{t,to) = H{t)U{t,to). (58) 

Further taking the initial condition 

U(t ,t )=I, (59) 

the evolution operator can be condensed into a single integral equation: 

U(t, t ) = I - - / H(t')U(t\ t )dt'. (60) 

When the operator H does not depend on time, Eq. ( 1601) can easily be integrated and at 
finally gives out 

U(t } t ) = e- iHit - to)/h . (61) 

2. Step by step time- evolution matrix transformation 

To unveil the behavior of concurrence at time t, we need to find the density matrix of 
the system at that moment, which can be obtained from 

p(t) = U(t)p(0)U\t). (62) 

Although Eq. (160]) gives a beautiful expression for the evolution operator, in reality U is 
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h(t) 



FIG. 22: The external magnetic field in a step function form h(t) = a @ t < to, h(t) = b @ t > to. 

hard to be obtained because of the integration involved. In order to overcome this obstacle, 
let us first consider the simplest time-dependent magnetic field: a step function of the form 

(Fig. 122) 

h(t) = a + (b - a)9(t - t ) , (63) 

where 9(t — to) is the usual mathematical step function defined by 

(64) 




At t and before, the system is time-independent since H a = H(t < to) — — J2<i ,> °f a j ~ 
aEiaf. Therefore we are capable of evaluating its ground state and density matrix at 
to straightforwardly. For the interval to to t, the Hamiltonian Hb = H(t > to) = 
~J2<ij>< J i°~ : j ~~ bTiiaf does not depend on time either, so Eq. ( IBTl) enables us to write 
out 

U(t,t )=e- iH ^ t>t ^ t - t0 ^ r \ (65) 

and therefore 

p(t) = ^(t,to)p(t )t/ t (t,to). (66) 

Starting from here, it is not hard to think of breaking an arbitrary magnetic function into 
small time intervals, and treating every neighboring intervals as a step function. Comparing 
the two graphs in Fig. [23], the method has just turned a ski sliding into a mountain climbing. 
Assuming each time interval is At, setting H — 1 then 

UfatoMo) = U(t il t i „ 1 )U(t i _ u t l _ 2 )...U(t u t )\tjj }, (67) 

i 

U(ti,t ) = Y[exp[-iH{t k )At], (68) 

fc=i 
U(U,t ) = exp[-iH(U)At]U(U-i-to)' (69) 
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FIG. 23: Divide an arbitrary magnetic field function into small time intervals. Every time step is At. 
Treat the field within the interval as a constant. In all, we turn a smooth function into a collection of step 
functions, which makes the calculation of dynamics possible. 

Here we avoided integration, instead we have chain multiplications which can be easily 
realized as loops in computational calculations. This is a common numerical technique; 
desired precisions can be achieved via proper time step length adjustment. 

3. Step by step projection 



Step by step matrix transformation method successfully breaks down the integration, 
but still involves matrix exponential, which is numerically resource costing. We propose 
a projection method to accelerate the calculations. Let us look at the step magnetic field 
again Fig. |22] For H a , after enough long time, the system at zero temperature is in the 
ground state |0) with energy, say, e. We want to ask how will this state evolves after the 
magnetic field is turned to the value b? Assuming the new Hamiltonian Hb has N eigenpairs 
Ei and \i/ji). The original state |0) can be expanded in the basis {|V'i)} : 



CilV'i) + c 2 \ip2) + ••• + c N \ij N ), 



(70) 



where 

Ci = (ipi\4>). 

When H is independent of time between t and to then we can write 

U(t, t ) |VW = e-^>^-">)/%, t0 > = e-^-^/%, i0 >, 



(71) 



(72) 
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Now the exponent in the evolution operator is a number no longer a matrix. The ground 
state will evolve with time as 

|0(t)) = cil^e-"^-* ) + c 2 \^ 2 )e- lE ^- to) + ... + c JV |V'iv)e-^ ( *-* o) 

= ^Cil^>e _iB,(t_t0) - (73) 

i=i 

and the pure state density matrix becomes 

p(t) = |0(t))(0(t)|. (74) 

Again any complicated function can be treated as a collection of step functions. When 
the state evolves to the next step just repeat the procedure to get the following results. Our 
test shows, for the same magnetic field both methods give the same results, but projection 
is much (about 20 times faster) than matrix transformation. This is a great advantage when 
the system size increases. But this is not the end of the problem. The summation is over 
all the eigenstates. Extending one layer out to 19 sites, fully diagonalizing the 2 19 by 2 19 
Hamiltonian and summing over all of them in every time step is breath taking. 

4- Dynamics of the spin system in a time- dependent magnetic field 

We consider the dynamics of entanglement in a two-dimensional spin system, where spins 
are coupled through an exchange interaction and subject to an external time-dependent mag- 
netic field. Four forms of time-dependent magnetic field are considered: step, exponential, 
hyperbolic and periodic. 

h{t) = a + (b - a)9(t - t ) , (75) 

.a t < tn 

h(t) ={ - (76) 

b + (a - 6)e- wt t > t 



(b-a) 



[tanh(wt) + 1] + a t > t 



, a t < tn 

h(t) = { " (78) 

a — a sm(ut + <fi) t > t 

We show in the following figures that the system entanglement behaves in an ergodic way 

in contrary to the one-dimensional Ising system. The system shows great controllability 
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under all forms of external magnetic field except the step function one which creates rapidly 
oscillating entanglement. This controllability is shown to be breakable as the different mag- 
netic field parameters increase. Also it will be shown that the mixing of even a few excited 
states by small thermal fluctuation is devastating to the entanglement of the ground state of 
the system. These can be explained by the Fermi's golden rule and adiabatic approximation. 
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FIG. 24: Dynamics of C(l,4) (solid blue line) in the 7-site system when the step magnetic field is changed 
from a = 0.5 to b — 1.5 (before the "critical point" h — 1.64), from a = 2 to b = 3 (after) and from a = 0.5 
to b = 3 (big step cross the "critical point"), where time t is in the unit of J -1 , the dashed red line stand 
for the concurrence corresponding to a constant magnetic field h = a, the straight solid yellow line for h = b 
and the dot-dashed green line for the average value of the oscillating concurrence. 



C. Tuning entanglement and ergodicity of spin systems using impurities and 
anisotropy 



There has been great interest in studying the differ ent sourc es of errors in quantum com- 
puting and their effect on quantum gate operations |l3ll Il32| . Different approaches have 
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FIG. 25: Dynamics of the concurrences C(l,2) (solid red (upper) line) and C(l,4) (solid blue (lower) line) 
in applied exponential magnetic fields of various frequencies u = 0.1, 0.5 & 1, with strength a = 1, b = 2, 
where time t is in the unit of J -1 . The straight dotted red (upper) lines are concurrences G(l,2) under 
constant magnetic field a = 1 and 6 = 2, so are the blue (lower) lines for C(l,4). 



been proposed for protecting quantum syste ms du ring the computational implementation of 
algorithms such as quantum error correction 123| and decoherence-free subspace 35| . Nev- 
ertheless, realizing a practical protection against the different types of induced decoherence 
is still a hard task. Therefore, studying the effect of naturally existing sources of errors 
such as impurities and lack of isotropy in coupling between the quantum systems imple- 
menting the quantum computing algorithms is a must. Furthermore, considerable efforts 
should be devoted to utilizing such sources to tune the entanglement rather than eliminat- 
ing them. The effect of impurities and anisotropy of coupling between neighbor spins in 
a one dimensional spin system has been investigated 46]. It was demonstrated that the 
entanglement can be tuned in a class of one- dimensional systems by varying the anisotropy 
of the coupling parameter as well as by introducing impurities into the spin system. For a 



18 




0.2r 



(0=0.1 




20 40 60 80 100 

t 



(a) 



(b) 



0.2 r 



0.2 r 



co=0.5 



(0=1 





40 60 80 

t 



100 



20 



40 60 80 

t 



100 



(c) 



(d) 



FIG. 26: Dynamics of the concurrences C(l,2) (solid red (upper) line) and C(l,4) (solid blue (lower) line) 
in applied hyperbolic magnetic fields of various frequencies u — 0.1, 0.5 & 1, with strength a = 1, b = 2, 
where time t is in the unit of J -1 . The straight dotted red (upper) lines are concurrences G(l,2) under 
constant magnetic field a = 1 and 6 = 2, so are the blue (lower) lines for C(l,4). 

physical quantity to be eligible for an equilibrium statistical mechanical description it has to 
be ergodic, which means that its time average coincides with its ensemble average. To test 
ergodicity for a physical quantity one has to compare the time evolution of its physical state 
to the corresponding equilibrium state. There has been an intensive efforts to investigate 
ergodicity in one-dimensional spin chains where it was demonstrated that the entanglement, 
magnetization, spin-spin correlation functions are non-ergodic in Isi ng and XY spi n chains 



for finite number of spins as well as at the thermodynamic limit 90 



111 



114 



133] 



In this part, we consider the entanglement in a two-dimensional XY triangular spin sys- 
tem, where the nearest neighbor spins are coupled through an exchange interaction J and 
subject to an external magnetic field h. We consider the system at different degrees of 
anisotropy to test its effect on the system entanglement and dynamics. The number of spins 
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FIG. 27: Dynamics of the concurrences C(l,2) (red (upper) line) and C(l,4) (blue (lower) line) in applied 
sine magnetic fields of various frequencies and field strength w = 0.1&a=l,w = 0.5&:a=l and u — 0.5 
& a = 5, where time t is in the unit of J -1 . 

in the system is 7, where all of them are identical except one (or two) which are considered 
impurities. The Hamiltonian for such a system is given by 

(1 + 7) y T _ x _ x (1-7) 



H 



Yl J ^ a i a l 



X)jyotoJ-M*)E°? 



(79) 



<»,j> 



where < i,j > is a pair of nearest-neighbors sites on the lattice, Jij = J for all sites 
except the sites nearest to an impurity site. For a single impurity, the coupling between 
the impurity and its neighbors Jij = J' = (a + 1) J, where a measures the strength of 
the impurity. For double impurities Jy = J' = («i + 1) J is the coupling between the two 
impurities and Jy = J" = (ct2 + 1) J is the coupling between any one of the two impurities 
and its neighbors while the coupling is just J between the rest of the spins. For this model 
it is convenient to set J = 1. Exactly solving Schrodinger equation of the Hamiltonian (1791 . 
yielding the system energy eigenvalues Ei and eigenfunctions ipi . The density matrix of the 
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system is defined by 

P=W)(M, (80) 

where |^o) is the ground state energy of the entire spin system. We confine our interest to 



the entanglement between two spins, at any sites i and j [72] . All the information needed in 
this case, at any moment t, is contained in the reduced density matrix Pij(t) which can be 
obtained from the entire system density matrix by integrating out all the spins states except 
% and j. We adopt the entanglement of form ation (or equivalently the concurrence C), as 



a well known measure of entanglement 109]. The dynamics of entanglement is eval uate d 



using the the step-by-step time-evolution projection technique introduced previously [134 ]. 



1. Single impurity 

We start by considering the effect of a single impurity located at the border site 1. The 
concurrence between the impurity site 1 and site 2, C(l,2), versus the parameter A for 
the three different models, Ising (7 = 1), partially anisotropic (7 = 0.5) and isotropic XY 
(7 = 0) at different impurity strengths (a = —0.5, 0, 0.5, 1) is in fig. [28j Firstly, the impurity 
parameter a is set to zero. For the corresponding Ising model, the concurrence C(l,2), in 
fig. [281(a). demonstrates the usual phase transition behavior where it starts at zero value and 
increases gradually as A increases reaching a maximum at A ~ 2 then decays as A increases 
further. As the degree of anisotropy decreases the behavior of the entanglement changes, 
where it starts with a finite value at A = and then shows a step profile for the small values 
of A. For the partially anisotropic case, the step profile is smooth and the entanglement 
mimics the Ising case as A increases but with smaller magnitude. The entanglement of 
isotropic XY system shows a sharp step behavior then suddenly vanishes before reaching 
A = 2. Interestingly, the entanglement behavior of the two-dimensional spin system at the 
different degrees of anisotropy mimics the behavior of the one-dimensional spin system at 
the same degrees of anisotropy at the extreme values of the parameter A. Comparing the 
entanglement behavior in the two-dimensional Ising spin system with the one-dimensional 
system, one can see a great resemblance except that the critical value becomes h/Jm2in 
the two dimensional case as shown in fig. [28j On the other hand, for the partially anisotropic 
and isotropic XY systems, the entanglements of the two-dimensional and one-dimensional 
system agrees at the extreme values of A where it vanishes for h » J and reaches a finite 
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FIG. 28: The concurrence (7(1, 2) versus the parameter A with a single impurity at the border site 1 with 
different impurity coupling strengths a — —0.5,0,0.5,1 for different degrees of anisotropy 7 = 1,0.5,0 as 
shown in the subfigures. The legend for all sub-figures is as shown in subfigure (a). 

value for h « J . The former case corresponds to an alignment of the spins in the in- 
direction, paramagnetic state, while the latter case corresponds to alignment in the x and 
y-directions which a ferromagnetic state. The effect of a weak impurity ( J' < J), a = —0.5, 
is shown in fig. 1281( b) where the entanglement behavior is the same as before except that 
the entanglement magnitude is reduced compared with the pure case. On the other hand, 
considering the effect of a strong impurity (J 7 > J), where a = 0.5 and 1, as shown in 
fig- [281( c) and fig. l28T d) respectively, one can see that the entanglement profile for 7 = 1 and 
0.5 have the same overall behavior as in the pure and weak impurity cases except that the 
entanglement magnitude becomes higher as the impurity gets stronger and the peaks shift 
toward higher A values. Nevertheless, the isotropic XY system behaves differently from the 
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previous cases where it starts to increase first in a step profile before suddenly dropping to 
zero again, which will be explained latter. 

To explore the effect of the impurity location we investigate the case of a single impurity 
spin located at site 4, instead of site 2, where we plot the concurrences C(l, 2) in figs. [2HJ 
Interestingly, while changing the impurity location has almost no effect on the behavior of 
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FIG. 29: The concurrence C(l,2) versus the parameter A with a single impurity at the central site 4 with 
different impurity coupling strengths a = —0.5,0,0.5,1 for different degrees of anisotropy 7 = 1,0.5,0 as 
shown in the subfigures. The legend for all subfigures is as shown in subfigure (a). 



the entanglement C(l,2) of the partially anisotropic and isotropic XY systems, it has a 
great impact on that of the Ising system where the peak value of the entanglement increases 
significantly in the weak impurity case and decreases as the impurity gets stronger as shown 
in fig. El 

Now we turn to the dynamics of the two dimensional spin system under the effect of a 
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FIG. 30: Dynamics of the concurrences C(l, 2), C(l, 4), C(2, 4) with a single impurity at the border site 1 
with different impurity coupling strengths a = —0.5, 0, 1, 2 for the two dimensional Ising lattice (7 = 1) under 
the effect of an exponential magnetic held with parameters values a=l, b=3.5 and u) = 0.1. The straight 
dashed lines represent the equilibrium concurrences corresponding to constant magnetic held h = 3.5. The 
legend for all subhgures is as shown in subhgure (a). 



single impurity and different degrees of anisotropy. We investigate the dynamical reaction 
of the system to an applied time-dependent magnetic field with exponential form h(t) = 
b + (a — b)e~ wt for t > and h(t) = a for t < 0. We start by considering the Ising 
system, 7 = 1 with a single impurity at the border site 1, which is explored in fig. [30] For 
the pure case, a = shown in fig. l3"UT a). the results confi rms the ergodic behavior of the 
system that was demonstrated in our previous work 134], where the asymptotic value of 
the entanglement coincide with the equilibrium state value at h(t) = b. As can be noticed 
from figs. 130T b). l30T c) and l30T d) neither the weak nor strong impurities have effect on the 
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ergodicity of the Ising system. Nevertheless, there is a clear effect on the asymptotic value 
of entanglements C(l, 2) and C(l, 4) but not on C(2, 4) which relates two regular sites. The 
weak impurity, a = —0.5 reduces the asymptotic value of C(l,2) and C(l,4) while the 
strong impurities, a = 1, 2 raise it compared to the pure case. 
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FIG. 31: Dynamics of the concurrences C(1,2),C(1,4),C(2,4) with a single impurity at the border site 
1 with different impurity coupling strengths a = —0.5,0,1,2 for the two dimensional partially anisotropic 
lattice (7 = 0.5) under the effect of an exponential magnetic field with parameters values a=l, b=3.5 
and u> = 0.1. The straight dashed lines represent the equilibrium concurrences corresponding to constant 
magnetic field h = 3.5. The legend for all subfigures is as shown in subfigure (b). 



The dynamics of the partially anisotropic XY system under the effect exponential mag- 
netic field with parameters a = 1, b = 3.5 and u = 0.1, is explored in fig. |3U It is remarkable 
to see that while for both the pure and weak impurity cases, a = and —0.5 , the system is 
nonergodic as shown in figs. l3Tl a) and l3~lTb). and it is ergodic in the strong impurity cases 
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a 



1 and 2 as illustrated in figs. l3lT c) and I3l7 d). 



2. Double impurities 







FIG. 32: The concurrence C(l,2), C(l,4), C(4, 5) versus the impurity coupling strengths ot\ and «2 with 
double impurities at sites 1 and 2 for the two dimensional Ising lattice (7 = 1) in an external magnetic field 
h=2. 

In this subsection we study the effect of double impurity, where we start with two located 
at the border sites 1 and 2. We set the coupling strength between the two impurities as 
J' = (1 + ai)J, between any one of the impurities and its regular nearest neighbors as 
J" = (1 + a-i) J and between the rest of the nearest neighbor sites on the lattice as J . The 
effect of the impurities strength on the concurrence between different pairs of sites for the 
Ising lattice is shown in fig. |32j In fig. 1327 a) we consider the entanglement between the two 
impurity sites 1 and 2 under a constant external magnetic field h = 2. The concurrence 
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C(l, 2) takes a large value when the impurity strengths ot,\, controlling the coupling between 
the impurity sites, is large and when a 2 , controlling coupling between impurities and their 
nearest neighbors, is weak. As a\ decreases and a 2 increases, C(l, 2) decreases monotonically 
until it vanishes. As one can conclude, a,\ is more effective than a 2 in controlling the 
entanglement in this case. On the other hand, the entanglement between the impurity 
site 1 and the regular central site 4 is illustrated in fig. 1327 b) which behaves completely 
different from C(l,2). The concurrence C(l,4) is mainly controlled by the impurity strength 
a 2 where it starts with a very small value when the impurity is very weak and increases 
monotonically until it reaches a maximum value at a 2 = 0, i.e. with no impurity, and 
decays again as the impurity strength increases. The effect of ot\ in that case is less significant 
and makes the concurrence slowly decreases as a,\ increases which is expected since as the 
coupling between the two border sites 1 and 2 increases the entanglement between 1 and 4 
decreases. It is important to note that in general C(l,2) is much larger than C(l,4) since 
the border entanglement is always higher than the central one as the entanglement is shared 
by many sites. The entanglement between two regular sites is shown in fig. l32T c) where the 
concurrence C(4,5) is depicted against a\ and a 2 , the entanglement decays gradually as a 2 
increases while a,\ has a very small effect on the entanglement, which slightly decreases as 
«i increases as shown. Interestingly, the behavior of the energy gap between the ground 
state and the first excited state of the Ising system AE versus the impurity strengths a,\ 
and a 2 , which is explored in fig. l32T d) has a strong resemblance to that of the concurrence 
C(4, 5) except that the decay of AE against a 2 is more rapid. The partially anisotropic 
system, 7 = 0.5, with double impurity at sites 1 and 2 and under the effect of the external 
magnetic field h — 2 is explored in fig. [33j As one can see, the overall behavior specially 
at the border values of the impurity strengths is the same as observed in the Ising case 
except that the concurrences suffer a local minimum within a small range of the impurity 
strength a 2 between and 1 while corresponding to the whole ot\ range. The change of the 
entanglement around this local minimum takes a step-like profile which is very clear in the 
case of the concurrence C(l,4) shown in fig. 133T b). Remarkably, the local minima in the 
plotted concurrences coincide with the line of vanishing energy gap as shown in fig. l33T d). 
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FIG. 33: The concurrence (7(1,2), (7(1,4), (7(4,5) versus the impurity coupling strengths ot\ and a.2 with 
double impurities at sites 1 and 2 for the two dimensional partially anisotropic lattice (7 = 0.5) in an 
external magnetic field h=2. 

3. Entanglement and quantum phase transition 

Critical quantum behavior in a many body system happens either when an actual crossing 
takes place between the excited state and the ground state or a limiting avoided level- crossing 
between them exists, i.e. an energy gab between the two states that vanishes in the infinite 

|. When a many body system crosses a critical 



system size limit at the critical point 
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point, significant changes in both its wave function and ground state energy takes place, 
which are manifested in the behavior of the entanglement function. The entanglement in 
one dimensional infinite spin systems, Ising and XY, was shown to demonstrate scaling 
behavior in the vicinity of critical points [72]. The change in the entanglement across the 
critical point was quantified by considering the derivative of the concurrence with respect 
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FIG. 34: (a) The concurrence Cu versus A; (b) the first derivative of the concurrence Cu with respect 
to A versus A; (c) the energy gap between the ground state and first excited state versus A; (d) the first 
derivative (in units of J) and second derivative (in units of J 2 ) of the energy gap with respect A versus A 
for the pure Ising system (7 = 1 and a = 0). 



to the parameter A. This derivative was explored versus A for different system sizes and 
although it didn't show divergence for finite system sizes, it showed clear anomalies which 
developed to a singularity at the thermodynamic limit. The ground state of the Heisenberg 
spin model is known to have a double degeneracy for an odd number of spins which is 
never achieved unless the thermodynamic limit is reached 68(. Particularly, the Ising ID 
spin chain in an external transverse magnetic field has doubly degenerate ground state in 
a ferromagnetic phase that is gapped from the excitation spectrum by 2J(1 — h/J), which 
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FIG. 35: (a) The concurrence Cu versus A; (b) the first derivative of the concurrence C14 with respect 
to A versus A; (c) the energy gap between the ground state and first excited state versus A; (d) the first 

second derivative (in units of J 2 ) of the energy gap with respect A versus A 
= 0.5 and a = 0). Notice that the first derivative of energy gap 
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is enlarged 10 times its actual scale for clearness. 



is removed at the critical point and the system becomes in a paramagnetic phase. Now let 
us first consider our two-dimensional finite size Ising spin system. The concurrence C14 and 
its first derivative are depicted versus A in figs. 1347 a) and (b) respectively. As one can see, 
the derivative of the concurrence shows strong tendency of being singular at A c = 1.64. The 
characteristics of the energy gap between the ground state and the first excited state as a 
function of A are explored in fig. l34T c). The system shows strict double degeneracy, zero 
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FIG. 36: (a) The concurrence Cu versus A; (b) the first derivative of the concurrence C14 with respect 
to A versus A; (c) the energy gap between the ground state and first excited state versus A; (d) the first 

second derivative (in units of J 2 ) of the energy gap with respect A versus A 
Notice that the first derivative of energy gap is enlarged 



derivative (in units of J) and 

for the pure isotropic system (7 = and a — 0). 

20 times its actual scale for clearness. 



energy gap, only at A = i.e. at zero magnetic field, but once the magnetic field is on the 
degeneracy is lifted and an extremely small energy gap develops, which increase very slowly 
for small magnetic field values but increases abruptly at certain A value. It is important to 
emphasis here that at A = 0, regardless of which one of the double ground states is selected 
for evaluating the entanglement, the same value is obtained. The critical point of a phase 
transition should be characterized by a singularity in the ground state energy, and an abrupt 
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change in the energy gap of the system as a function of the system parameter as it crosses the 
critical point. To better understand the behavior of the energy gap across the prospective 
critical point and identify it, we plot the first and second derivatives of the energy gap as a 
function of A in fig. IMT d). Interestingly, the first derivative dAE/dX which represents the 
rate of change of the energy gap as a function of A starts with a zero value at A = and then 
increase very slowly before it shows a great rate of change and finally reaches a saturation 
value. This behavior is best represented by the second derivative d 2 AE/d\ 2 , which shows 
strong tendency of being singular at A c = 1.8, which indicates the highest rate of change the 
energy gap as a function of A. The reason for the small discrepancy between the two values 
of the A c extracted from the dC/dX plot and the one of d 2 AE/dX 2 is that the concurrence 
C14 is only between two sites and does not represent the whole system in contrary to the 
energy gap. One can conclude that the rate of change of the energy gap as a function of 
the system parameter, A in our case, should be maximum across the critical point. Turning 
to the case of the partially anisotropic spin system, 7 = 0.5, presented in fig. [35j one can 
notice from fig. l35Ta) that the concurrence shows few sharp changes, which is reflected in 
the energy gap plot as an equal number of minima as shown in fig. 135T b). Nevertheless, 
again there is only one strict double degeneracy at A = while the other three energy gap 
minima are non-zero and in the order of 10 -5 . It is interesting to notice that the anomalies 
in both dC/dX and d 2 AE/dX 2 are much stronger and sharper compared with the Ising case 
as shown in figs. l35Tc) and (d). Finally the isotropic system which is depicted in fig. EH 
shows even sharper energy gap changes as a result of the sharp changes in the concurrence 
and the anomalies in the derivatives dC/dX and d 2 AE/dX 2 are even much stronger than the 
previous two cases. 
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